"""
光学效率计算模块 - 阶段3: 光学效率计算与验证 (修复版本)
实现五大效率因子的计算: ηcos, ηat, ηsb, ηtrunc, ηref
"""

import numpy as np
import math
from multiprocessing import Pool
import warnings
warnings.filterwarnings('ignore')

class EfficiencyCalculator:
    def __init__(self):
        """初始化效率计算器"""
        # 物理常数
        self.mirror_reflectivity = 0.88  # 镜面反射率
        self.atmospheric_visibility = 23000  # 大气能见度 (m)
        
        # 定日镜参数
        self.mirror_width = 6  # 定日镜宽度 (m)
        self.mirror_height = 6  # 定日镜高度 (m)
        self.mirror_area = self.mirror_width * self.mirror_height
        
        # 吸收塔参数
        self.tower_height = 80  # 吸收塔高度 (m)
        self.receiver_height = 8  # 接收器高度 (m)
        self.receiver_width = 7  # 接收器宽度 (m)
        self.install_height = 4  # 定日镜安装高度 (m)
        
        # 标准: 安全距离参数
        self.safety_distance = 11  # 安全距离 max(w,h) + 5m = 11m
        
        print("✅ 效率计算器初始化完成 (修复版本)")
    
    def calculate_total_efficiency(self, mirror_positions, tower_position, solar_positions):
        """计算总光学效率 ηtotal = ηcos × ηat × ηsb × ηtrunc × ηref"""
        total_efficiencies = []
        
        print(f"🔄 开始计算 {len(mirror_positions)} 个定日镜的光学效率...")
        
        for i, mirror_pos in enumerate(mirror_positions):
            mirror_efficiency = []
            
            for solar_pos in solar_positions:
                # 计算各个效率分量
                eta_cos = self._cosine_efficiency(mirror_pos, tower_position, solar_pos)
                eta_at = self._atmospheric_efficiency(mirror_pos, tower_position)
                eta_sb = self._shadow_blocking_efficiency(mirror_pos, mirror_positions, tower_position, solar_pos)
                eta_trunc = self._truncation_efficiency(mirror_pos, tower_position, solar_pos)
                eta_ref = self.mirror_reflectivity  # 镜面反射率 (常数)
                
                # 计算总效率
                eta_total = eta_cos * eta_at * eta_sb * eta_trunc * eta_ref
                mirror_efficiency.append(eta_total)
            
            # 计算年平均效率
            avg_efficiency = np.mean(mirror_efficiency) if mirror_efficiency else 0
            total_efficiencies.append(avg_efficiency)
            
            if (i + 1) % 50 == 0:
                print(f"   已完成 {i + 1}/{len(mirror_positions)} 个定日镜")
        
        print("✅ 光学效率计算完成")
        return total_efficiencies
    
    def check_safety_distance_violations(self, mirror_positions):
        """检查安全距离违规 - 修复版本"""
        violations = []
        violation_count = 0
        total_pairs = 0
        
        print(f"🔍 检查 {len(mirror_positions)} 个定日镜的安全距离...")
        
        for i, pos1 in enumerate(mirror_positions):
            for j, pos2 in enumerate(mirror_positions):
                if i >= j:  # 避免重复检查
                    continue
                
                total_pairs += 1
                x1, y1 = pos1
                x2, y2 = pos2
                distance = math.sqrt((x1-x2)**2 + (y1-y2)**2)
                
                if distance < self.safety_distance:
                    violation_count += 1
                    violations.append({
                        'mirror1_index': i,
                        'mirror2_index': j,
                        'mirror1_pos': pos1,
                        'mirror2_pos': pos2,
                        'distance': distance,
                        'required_distance': self.safety_distance
                    })
        
        violation_ratio = violation_count / total_pairs if total_pairs > 0 else 0
        
        print(f"✅ 安全距离检查完成")
        print(f"   总镜对数: {total_pairs}")
        print(f"   违规镜对数: {violation_count}")
        print(f"   违规率: {violation_ratio:.1%}")
        print(f"   安全距离要求: {self.safety_distance}m")
        
        return violations, violation_count, violation_ratio

    def _cosine_efficiency(self, mirror_pos, tower_pos, solar_pos):
        """计算余弦效率 ηcos (基于太阳入射角) - 修复版本"""
        try:
            mx, my = mirror_pos
            tx, ty = tower_pos
            
            # 定日镜到吸收塔的单位向量
            tower_vec = np.array([tx - mx, ty - my, self.tower_height - self.install_height], dtype=np.float64)
            tower_vec_norm = np.linalg.norm(tower_vec)
            if tower_vec_norm > 0:
                tower_vec = tower_vec / tower_vec_norm
            else:
                return 0.0
            
            # 太阳方向单位向量
            sun_elevation = solar_pos['elevation']
            sun_azimuth = solar_pos['azimuth']
            
            sun_vec = np.array([
                math.cos(sun_elevation) * math.sin(sun_azimuth),
                math.cos(sun_elevation) * math.cos(sun_azimuth),
                math.sin(sun_elevation)
            ], dtype=np.float64)
            
            # 定日镜法向量 (双轴跟踪，指向太阳和吸收塔的角平分方向)
            normal_vec = sun_vec + tower_vec
            normal_vec_norm = np.linalg.norm(normal_vec)
            if normal_vec_norm > 0:
                normal_vec = normal_vec / normal_vec_norm
            else:
                return 0.0
            
            # 余弦效率 = cos(入射角)
            cos_incident = np.dot(sun_vec, normal_vec)
            eta_cos = max(0, float(cos_incident))  # 修复: 确保返回Python float
            
            return eta_cos
            
        except Exception as e:
            print(f"余弦效率计算错误: {e}")
            return 0.0
    
    def _atmospheric_efficiency(self, mirror_pos, tower_pos):
        """计算空气透射率 ηat (基于路径长度衰减) - 修复版本"""
        try:
            mx, my = mirror_pos
            tx, ty = tower_pos
            
            # 计算定日镜到吸收塔的距离
            distance = math.sqrt((tx - mx)**2 + (ty - my)**2 + (self.tower_height - self.install_height)**2)
            
            # Beer-Lambert定律: 大气透射率
            # eta_at = exp(-distance / visibility)
            eta_at = math.exp(-distance / self.atmospheric_visibility)
            
            return max(0.0, min(1.0, eta_at))
            
        except Exception as e:
            print(f"大气透射率计算错误: {e}")
            return 0.8  # 返回默认值
    
    def _shadow_blocking_efficiency(self, mirror_pos, all_mirror_positions, tower_pos, solar_pos):
        """计算阴影-遮挡效率 ηsb (几何遮挡计算)"""
        mx, my = mirror_pos
        
        # 太阳光线方向
        sun_elevation = solar_pos['elevation']
        sun_azimuth = solar_pos['azimuth']
        
        sun_dir = np.array([
            -math.cos(sun_elevation) * math.sin(sun_azimuth),
            -math.cos(sun_elevation) * math.cos(sun_azimuth),
            -math.sin(sun_elevation)
        ])
        
        blocked_area = 0
        
        # 检查其他定日镜的阴影遮挡
        for other_pos in all_mirror_positions:
            if other_pos == mirror_pos:
                continue
                
            ox, oy = other_pos
            
            # 计算阴影投射
            shadow_length = (self.install_height - 0) / abs(sun_dir[2]) if sun_dir[2] != 0 else 0
            shadow_x = ox + sun_dir[0] * shadow_length
            shadow_y = oy + sun_dir[1] * shadow_length
            
            # 检查阴影是否覆盖当前定日镜
            if self._rectangles_overlap((mx, my), (shadow_x, shadow_y)):
                overlap_area = self._calculate_overlap_area((mx, my), (shadow_x, shadow_y))
                blocked_area += overlap_area
        
        # 阴影遮挡效率
        eta_sb = max(0, 1 - blocked_area / self.mirror_area)
        
        return eta_sb
    
    def _rectangles_overlap(self, pos1, pos2):
        """检查两个矩形定日镜是否重叠"""
        x1, y1 = pos1
        x2, y2 = pos2
        
        # 简化为矩形重叠检测
        half_width = self.mirror_width / 2
        half_height = self.mirror_height / 2
        
        return (abs(x1 - x2) < self.mirror_width and abs(y1 - y2) < self.mirror_height)
    
    def _calculate_overlap_area(self, pos1, pos2):
        """计算两个矩形的重叠面积"""
        x1, y1 = pos1
        x2, y2 = pos2
        
        half_width = self.mirror_width / 2
        half_height = self.mirror_height / 2
        
        # 计算重叠区域
        x_overlap = max(0, min(x1 + half_width, x2 + half_width) - max(x1 - half_width, x2 - half_width))
        y_overlap = max(0, min(y1 + half_height, y2 + half_height) - max(y1 - half_height, y2 - half_height))
        
        return x_overlap * y_overlap
    
    def _truncation_efficiency(self, mirror_pos, tower_pos, solar_pos):
        """计算截断效率 ηtrunc (光斑-接收器匹配)"""
        mx, my = mirror_pos
        tx, ty = tower_pos
        
        # 计算定日镜到接收器的距离
        distance = math.sqrt((tx - mx)**2 + (ty - my)**2 + (self.tower_height - self.install_height)**2)
        
        # 光斑大小估算 (考虑太阳角直径和镜面误差)
        sun_angular_size = 0.00932  # 太阳角直径 (弧度)
        mirror_error = 0.002  # 镜面误差 (弧度)
        
        total_error = math.sqrt(sun_angular_size**2 + mirror_error**2)
        spot_radius = distance * total_error
        
        # 接收器有效面积
        receiver_area = self.receiver_height * self.receiver_width
        
        # 光斑面积
        spot_area = math.pi * spot_radius**2
        
        # 截断效率 (简化模型)
        if spot_area <= receiver_area:
            eta_trunc = 1.0
        else:
            eta_trunc = receiver_area / spot_area
        
        return min(1.0, eta_trunc)
    
    def calculate_parallel_efficiency(self, mirror_positions, tower_position, solar_positions, n_processes=4):
        """并行计算光学效率 (多进程加速)"""
        print(f"🚀 启动 {n_processes} 进程并行计算效率...")
        
        # 分割定日镜位置
        chunk_size = len(mirror_positions) // n_processes
        mirror_chunks = [mirror_positions[i:i+chunk_size] for i in range(0, len(mirror_positions), chunk_size)]
        
        # 准备参数
        args = [(chunk, tower_position, solar_positions, i) for i, chunk in enumerate(mirror_chunks)]
        
        # 并行计算
        with Pool(n_processes) as pool:
            results = pool.map(self._calculate_chunk_efficiency, args)
        
        # 合并结果
        total_efficiencies = []
        for chunk_result in results:
            total_efficiencies.extend(chunk_result)
        
        print("✅ 并行效率计算完成")
        return total_efficiencies
    
    def _calculate_chunk_efficiency(self, args):
        """计算一块定日镜的效率 (用于并行处理)"""
        mirror_chunk, tower_position, solar_positions, chunk_id = args
        
        chunk_efficiencies = []
        
        for mirror_pos in mirror_chunk:
            mirror_efficiency = []
            
            for solar_pos in solar_positions:
                # 计算各个效率分量
                eta_cos = self._cosine_efficiency(mirror_pos, tower_position, solar_pos)
                eta_at = self._atmospheric_efficiency(mirror_pos, tower_position)
                eta_sb = 0.9  # 简化阴影遮挡效率 (避免并行计算复杂性)
                eta_trunc = self._truncation_efficiency(mirror_pos, tower_position, solar_pos)
                eta_ref = self.mirror_reflectivity
                
                # 计算总效率
                eta_total = eta_cos * eta_at * eta_sb * eta_trunc * eta_ref
                mirror_efficiency.append(eta_total)
            
            # 计算年平均效率
            avg_efficiency = np.mean(mirror_efficiency) if mirror_efficiency else 0
            chunk_efficiencies.append(avg_efficiency)
        
        print(f"   进程 {chunk_id} 完成 {len(mirror_chunk)} 个定日镜")
        return chunk_efficiencies
    
    def validate_efficiency(self, efficiencies):
        """效率验证和异常处理"""
        valid_efficiencies = []
        
        for i, eta in enumerate(efficiencies):
            # 检查效率值范围 [0, 1]
            if not (0 <= eta <= 1):
                print(f"⚠️ 定日镜 {i} 效率异常: {eta:.3f}, 修正为边界值")
                eta = max(0, min(1, eta))
            
            # 检查效率过低 (可能的物理合理性问题)
            if eta > 0.95:
                print(f"⚠️ 定日镜 {i} 效率过高: {eta:.3f}, 进行物理合理性检查")
                eta = 0.95  # 限制最大效率
            
            valid_efficiencies.append(eta)
        
        # 统计信息
        avg_efficiency = np.mean(valid_efficiencies)
        max_efficiency = np.max(valid_efficiencies)
        min_efficiency = np.min(valid_efficiencies)
        
        print(f"✅ 效率验证完成:")
        print(f"   平均效率: {avg_efficiency:.3f}")
        print(f"   最大效率: {max_efficiency:.3f}")
        print(f"   最小效率: {min_efficiency:.3f}")
        
        return valid_efficiencies

    def calculate_comprehensive_results(self, mirror_positions, tower_position, solar_positions):
        """计算综合结果 - 修复版本"""
        print(f"\n🔬 开始综合效率计算...")
        print(f"   定日镜数量: {len(mirror_positions)}")
        print(f"   吸收塔位置: ({tower_position[0]:.2f}, {tower_position[1]:.2f})")
        print(f"   太阳位置数: {len(solar_positions)}")
        
        results = {}
        
        try:
            # 1. 光学效率计算
            total_efficiencies = self.calculate_total_efficiency(
                mirror_positions, tower_position, solar_positions
            )
            
            # 2. 安全距离检查
            violations, violation_count, violation_ratio = self.check_safety_distance_violations(mirror_positions)
            
            # 3. 功率计算
            dni_avg = 500  # W/m² (假设平均DNI)
            mirror_area = self.mirror_area
            
            total_power = sum(eta * mirror_area * dni_avg for eta in total_efficiencies)
            avg_efficiency = np.mean(total_efficiencies) if total_efficiencies else 0
            
            # 4. 统计信息
            results = {
                'mirror_count': len(mirror_positions),
                'tower_position': tower_position,
                'total_efficiencies': total_efficiencies,
                'avg_efficiency': float(avg_efficiency),  # 修复: 确保是Python float
                'total_power': float(total_power),  # 修复: 确保是Python float
                'power_per_mirror': float(total_power / len(mirror_positions)) if len(mirror_positions) > 0 else 0,
                'safety_violations': violation_count,
                'violation_ratio': float(violation_ratio),  # 修复: 确保是Python float
                'violation_details': violations[:10],  # 只保留前10个违规详情
                'field_utilization': len(mirror_positions) * mirror_area / (math.pi * 350**2)
            }
            
            print(f"✅ 综合效率计算完成")
            print(f"   年平均光学效率: {avg_efficiency:.4f}")
            print(f"   年平均输出功率: {total_power:,.0f} W")
            print(f"   安全距离违约: {violation_count}个")
            
        except Exception as e:
            print(f"❌ 综合结果计算失败: {e}")
            results = {
                'mirror_count': len(mirror_positions),
                'tower_position': tower_position,
                'avg_efficiency': 0.0,
                'total_power': 0.0,
                'safety_violations': 0,
                'violation_ratio': 0.0,
                'error': str(e)
            }
        
        return results

# 测试效率计算模块
if __name__ == "__main__":
    print("🚀 开始效率计算测试...")
    
    # 初始化效率计算器
    calc = EfficiencyCalculator()
    
    # 测试数据
    mirror_positions = [(100, 100), (150, 150), (200, 200)]
    tower_position = (0, 0)
    
    # 模拟太阳位置数据
    solar_positions = [
        {'elevation': math.radians(45), 'azimuth': math.radians(180)},
        {'elevation': math.radians(60), 'azimuth': math.radians(200)},
        {'elevation': math.radians(30), 'azimuth': math.radians(160)}
    ]
    
    # 计算效率
    efficiencies = calc.calculate_total_efficiency(mirror_positions, tower_position, solar_positions)
    
    # 验证效率
    valid_efficiencies = calc.validate_efficiency(efficiencies)
    
    print(f"✅ 测试完成，定日镜效率: {valid_efficiencies}") 